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Abstract 

The onset of the Rayleigh- Taylor instability is studied a compressible Brownian Yukawa fluid 
mixture on the "molecular" length and time scales of the individual particles. As a model, a two- 
dimensional phase-separated symmetric binary mixture of colloidal particles of type A and B with 
a fluid-fluid interface separating an j4-rich phase from a S-rich phase is investigated by Brownian 
computer simulations when brought into non-equilibrium via a constant external driving field 
which acts differently on the different particles and perpendicular to the interface. Two different 
scenarios are observed which occur either for high or for low interfacial free energies as compared to 
the driving force. In the first scenario for high interfacial tension, the critical wavelength A c of the 
unstable interface modes is in good agreement with the classical Rayleigh- Taylor formula provided 
dynamically rescaled values for the interfacial tension are used. The wavelength A c increases 
with time representing a self-healing effect of the interface due to a local density increase near 
the interface. The Rayleigh- Taylor formula is confirmed even if A c is of the order of a molecular 
correlation length. In the second scenario for very large driving forces as compared to the interfacial 
line tensions, on the other hand, the particle penetrate easily the interface by the driving field and 
form microscopic lanes with a width different from the predictions of the classical Rayleigh- Taylor 
formula. The results are of relevance for phase-separating colloidal mixtures in a gravitational or 
electric field. 
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I. INTRODUCTION 



Interface instabilities in driven non-equilibrium systems are well-known from macroscopic 
hydrodynamics . Examples include the classic Rayleigh- Taylor instability |l|of a heavy 
liquid on top of another lighter liquid, the Mullins-Sekerka fingering instability 3, 3] i n dif- 
fusive systems and the Saffmann- Taylor instability js,El for compressed liquids of different 
viscosities. However, what is much less clear is the microscopic origin of these instabilities, 
i.e. a view which resolves the microscopic discrete particle trajectories causing the insta- 
bility. Recently there has been progress in simulating large systems with 10000-10000000 
of discrete particles forming an interface Q, Isl 0, I n particular, the microscopic origin 
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of the Rayleigh- Taylor instability was explored in huge molecular dynamics p, |9|, |10j and 
direct numerical |8] computer simulations. A natural question concerns the applicability of 
coarse-grained hydrodynamics towards microscopic spatial dimensions. The basic quantities 
entering in hydrodynamics auch as the viscosity or the surface tension are quantities which 
are only defined for large systems and exhibit finite-size corrections when applied to small 
inhomogeneities. 

Another approach to interfacial instabilities is via mesoscopic colloidal suspensions which 
bear the fascinating possibility to study the particle trajectories in real-space by video mi- 
croscopy in quasi- two-dimensional suspensions |ll( . Well-characterized colloidal suspensions 
serve as model systems for many questions of many-body systems including inter- 

faces jl^, [3, Q, 12 , [3 [3| • The strength of an external field is easily tunable for colloids 
in striking contrast to molecular liquids. The dynamics of the mesoscopic colloidal parti- 
cles which are embedded in a microscopic solvent is Brownian rather than Newtonian [2^ 
such that dynamical quantities are different. To the best of our knowledge, for compressible 
Brownian fluids nothing is known about the onset of the Rayleigh- Taylor instability on a 



length scale of interparticle distances. Clearly, it will be different from mo 



ecular dynam- 



One 



ics where heat will be generated and inertia effects can lead to turbulence 
promising investigation in the Newtonian case is for a suspension of heavy granular grains 
where Rayleigh- Taylor instabilities have recently been observed by Rehberg and coworkers 
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I but also in this work the individual particle trajectories were not resolved. Mix- 
tures of colloids and polymers with a real-space analysis of the colloids trajectories represent 
another valuable system to look at for interfacial instabilities in external driving fields jl^j ]. 
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In the present paper we study, by Brownian dynamics computer simulations, the particle- 
resolved onset of the Rayleigh- Taylor instability. In our model, an interface separating an 
A-rich from an A-poor fluid is exposed to an external field which acts differently on the 
different particles and is directed perpendicular to the interface. In order to keep the model 
simple |23j and to link to two-dimensional colloidal suspensions, we take two spatial dimen- 
sions and a symmetric equimolar mixture interacting via a Yukawa pair potential. Most 
of our characteristics, however, will carry over to three-dimensional systems, to asymmetric 
mixtures, and to different interparticle interactions. 

The external field is so large that it will induce local density inhomogeneities such that 
the fluid is compressible in contrast to the ordinary case treated within the Navier-Stokes 
equations jicl |. We observe two different scenarios which occur either for high or for low 
interfacial free energies as compared to a typical thermal energy per average particle distance. 
In the first scenario for high interfacial tension, an interfacial instability is observed which is 
driven by the competition of the cost in interfacial tension versus gain in potential energy, 
similar in spirit to the classical Rayleigh- Taylor instability. The classical threshold value 
for the wavelength of unstable interface perturbations is confirmed provided the line tension 
is dynamically rescaled in terms of the actual local density at the interface. The critical 
wavelength separating stable from unstable density undulations increases with time. This 
can be interpreted as a "self-healing effect" of the interface caused by a local density increase 
near the interface. In the second scenario for small interfacial tensions, on the other hand, the 



particle penetrate easily the interface by t 
to previous simulation studies (24 , Q, U. 



re driving field and form microscopic lanes similar 



27, 



In this case the classical threshold for the 



unstable wavelength is smaller than a molecular correlation length such that a breakdown 
of macroscale hydrodynamics is expected. Our results are of relevance for phase-separating 
colloidal mixtures in a gravitational or electric field where similar effects have recently been 
reported by Aarts and coworkers Q| and for settling granular grains j^ . 

The paper is organized as follows: In section InJ we define the model used and describe 
briefly our simulation technique. As a prerequisite, a microscopic calculation of the density 
profile and the line tension is presented in section IIHI We further review the classical 
Rayleigh- Taylor instability briefly in section HVl Results for the interface instability for high 
line tensions are presented in section while the case of small line tensions is described 
in section IVII Conclusion are given in section IVHI In particular, we discuss a possible 



3 



verification of our predictions in experiments. 



II. THE MODEL 



In our model 



we consider a symmetric binary colloidal mixture comprising N 



Na + Nb Brownian colloidal particles in d = 2 spatial dimensions. The particles are in 



an area S with a fixed total number density of p = ^. Half of them are particles of type 



s 

A, the other half is of type B with partial number density Pa = Pb = f- The colloidal 
suspension is held at fixed temperature T via the bath of microscopic solvent particles. The 
colloidal particles of species a and b (a,b G {A, B}) are interacting via an effective Yukawa 
pair potential 

Vab(r) TT exp [-«(r - a ah )/a\ 

= U (Tab , (1) 

k B T r 

Here r is the center-to-center separation, ksT is the thermal energy, Uq is a dimensionless 

amplitude, a is the particle diameter as a length scale and k is the inverse screening length. 

The set of diameters a a b is taken as 



&AA = &BB = O (2) 



GAB = 0"(1 + A) (3) 

We chose A > corresponding to positive non-additivity. This implies that the cross- 
interaction Vab(j) is more repulsive than Vaa(j) = Vbb (y ) > which drives phase separation 
into an A-rich and a S-rich phase. 

The dynamics of the colloidal particles is overdamped Brownian motion. The friction 
constant £ = 37cr]a (with t] denoting the shear viscosity of the solvent) is assumed to be the 
same for both A and B particles. The constant external force for the ith particle of species 
a, F^ a \ is acting differently on the both constituents of the binary mixture. It is = e z F 
for A particles and F\ = —e z F for B particles. The stochastic Langevin equations for the 
colloidal trajectories f^ a \t) with i — 1, N a (a e {A, B}) read as 

= E^T a) + + Fi a) + ^(t). (4) 

3=1 3=1 
3¥"i 
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where 



r ij H ~ ^ \i an d ^ i s the complementary index to a (b = A if a = B and b = B if 
a = A). The right-hand-side includes all forces acting onto the colloidal particles, namely the 
force resulting from inter-particle interactions, the external constant force, and the random 
forces .Fj describing the collisions of the solvent molecules with the ith colloidal particle. 



i? 6) =-V^)K6(rJf), (5) 



The latter are Gaussian random numbers with zero mean, = 0, and variance 



(F^) a {t)(F^ ] )p(t') = 2k B T£6 a p6 tJ 6(t - t'). (6) 

The subscripts a and (3 stand for the two Cartesian components. In the absence of an 
external field and for A = 0, the model reduces to a two-dimensional Brownian Yukawa fluid 
in equilibrium which has been extensively investigated as far as structural and dynamical 
equilibrium correlations and freezing transitions are concerned (29 . 



30, 



3 it ] . For a positive A 



and vanishing external drive, our system will lead to equilibrium fluid-fluid phase separation 



including a critical point which has 



3een studied in non-additive hard-core models by theory 



and simulation, see e.g. Ref. 



We solve the Langevin equations of motion by Brownian dynamics simulations (BD) 
using a finite time-step At and the technique of Ermak |4l|, |42|. The typical 
size of the time-step was 0.003tb, where t b = £cr 2 /Vo is a suitable Brownian timescale. As 
a reference, the field- free case F = is studied extensively first. This is an equilibrium 
situation. In the simulation set-up here, we put Na = 1000 A and Ng = 1000 B particles 
into a rectangular cell of lengths L in x-direction and D in ^-direction with = | such that 
the total colloidal number density is p = -Mr. Then the system spontaneously exhibits a 
fluid-fluid interface along the x-direction separating an A-rich from a 5-rich fluid provided 
the density is larger than the critical density. To have an single interface system, without 
disturbing effects due to external walls, periodic boundary conditions are used in the re- 
direction while antiperiodic boundary conditions are used in z-direction. With antiperiodic 
boundary conditions, the particle type is changed from A to B if an A particle is crossing a 
boundary in z-direction and vice versa. In a finite slab around this boundary, all interactions 
are set to be equal, i.e. the cross interaction is Vab = Vaa = Vbb- This avoids a second 
interface and minimizes finite-size effects. Similar boundary condition shave been employed 
to study interfaces of symmetric polymer mixtures by Miiller, Binder and coworkers |43j . 



We equillibrated 3 x 10 5 time steps corresponds to an equilibration time of 900tb and 
gathered statistics again for a further simulation time of 900r B . We also checked that our 
data do not suffer from finite size effects by choosing a larger system with = 4000 A 



interface to be parallel to the x-axis. All our investigations in nonequilibrium are with this 
larger system size. 

With an equilibrated interface as a starting configuration we suddenly turn on the external 
field which drives the particles against the interface. An interface instability was observed. 
300 different equilibrated starting configuration were then typically used in order to perform 
time-dependent averages. For our nonequilibrium simulations, the antiperiodic boundary 
conditions in z direction has no consequences once the particles are driven against the 
interface. 

III. DENSITY PROFILES AND EQUILIBRIUM FLUID-FLUID INTERFACIAL 
FREE ENERGY 

In equilibrium, a phase-separated binary fluid mixture is characterized by its ^-dependent 
partial density profiles Pa( z ) and Pb{z) as defined via 



where a G {A, B} and (. . .) denotes a canonical average. Furthermore the associated fluid- 
fluid interfacial free energy 7 is a key quantity. Since our model is two-dimensional, this 
interfacial free energy 7 corresponds to a line tension. In a symmetric binary mixture, the 
partial density profiles and the line tension 7 depend on the temperature and the total 
number density p. 

Obviously, due to the A-B symmetry of our model, the density profiles are symmetric, 
i.e. pa{z) = pb{—z), if the interface position is at z = 0. Computer simulation results of 
the density profiles for a strong positive nonadditivity A = 1.6 and various bulk densities 
at a fixed temperature are presented in Figure 1. For high densities one clearly identifies a 
depletion zone in the interface as generated by the large non-additivity. For densities closer to 
the critical one, however, the density profiles are getting flatter and particles interpenetrate 
mutually. For large bulk densities, a slight density oscillation shows up as typical for fluid 



and N B = 4000 B particles with f 



The antiperiodic boundary conditions force the 




(7) 
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interfaces when the Fisher- Widom line is exceeded (44). Clearly, for the parameters chosen, 
the mixture is almost completely phase-separated, i.e. the partial densities of the A particles 
is practically zero in the B-region and vice-versa. Again this is different very close to the 
critical point. 

We have also calculated the equilibrium line tension as a function of density for fixed 



temperature by using "exact" computer simulation techniques. The interfacial 



can be gained by integration of the anisotropy of the pressure tensor (45 
two spatial dimensions, this expression reads 
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ree energy 
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J [P N (z) - P T (z)} dz 
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(8) 



where Pn(z) = p zz (z) and Pt{z, 

tensor. This local pressure tensor has the following form 



PxM are the diagonal ^onents ofthe loca! pressure 



1 componei 
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53. 54] 



Pap (z) = (p A (z) + p B (z))k B TS a0 




(9) 



l<i<N A 
l<j<N B 



where the subscripts a and (3 stand for the two Cartesian components, r. 



-(ab) _ -(a) -(b) 



'J 



r 



zf* — zf^ and Q(z) is the Heaviside step function. The z-dependence of the anisotropy 



of the pressure tensor is plotted in Figure 1 as well. As can be deduced from Figure 1, its 
main weight is centered in the interface position around z = 0. The oscillations of the 
density field for high densities correlate to that of the pressure tensor. 

Simulation results for the density- dependent line tension 'j(p) for two different nonaddi- 
tivities A = 0.8 and A = 1.6 are presented in Figure 2. Most of our calculations were for 
parameter combinations well away from the critical point where the line tension is mainly 
governed by internal energy such that the anisotropy of the pressure tensor is significant 
resulting in a relative small statistical error for 7, see Figure 2a. Some further points are 
also for smaller densities closer to the critical point, see Figure 2b. By crudely extrapolat- 
ing the data one could estimate the critical density to be at p c o 2 w 0.08 for A = 0.8 and 
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p c a 2 ~ 0.04 for A = 1.6 as indicated by arrows in Figure 2b. This extrapolation was cross- 
checked by mapping our Yukawa system onto that of effective non-additive hard disks using 
the prescription of Barker-Henderson for the effective hard-core diameter j^. The critical 
point can then be read off from that of a non-additive symmetric hard disk binary mixture 
which has been studied in detail via theory and simulation by Giaquinta and coworkers 36j . 
We find good agreement with our extrapolation as compared to the mapping procedure. 
As can be further deduced from Figure 2, the line tension "f(p) is strongly increasing with 
density. This is expected as a density increase means a smaller spacing between A and B 
particles such that the energetic non-additivity is getting more pronounced. At fixed density 
and temperature, the line tension increases with A which is clearly due to the fact that a 
larger nonadditivity leads to a larger energy cost of different particle species meeting at the 
interface. 

IV. CLASSICAL RAYLEIGH-TAYLOR INSTABILITY 

The classical Rayleigh- Taylor instability is obtained for a heavy incompressible liquid on 



top of a lighter incompressible liquid 



i. 



|56j. A small harmonic interfacial undulation with 
a wave length A yields a favorable decrease in potential energy but at the same time a free 
energy penalty due to the increasing arc length of the interface which costs line tension. If 
the wave length is larger than a critical one A > A c , however, the penalty is smaller than 
the gain such that an unstable mode whose amplitude grows in time is present. The critical 
wave length can be calculated in our two-dimensional situation as 



X c = 2 -^ = 2n L 7 . , (10) 

where 7 is the line tension as introduced in the previous section. Since the concept of line 
tension is a macroscopic one, Eq. fjlOj) is only justified as long as A c is much larger than any 
microscopic distance as e.g. the mean interparticle spacing a = 4=. This requires that the 
line tension has to be much larger than the driving force difference 



A c > a =^ 7 > 

In the opposite limit 



p{A) _ p(B) 



(11) 



7 « F(A)_P(B) (12) 
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it is expected that the classical concept of the Rayleigh- Taylor instability will break down. 
We shall explore this case in detail in section VI. Another question concerns the applicability 
of the Rayleigh- Taylor instability criterion when A c is of a similar order than the mean 
interparticle spacing a. As we shall show below, the criterion works remarkable well even if 
A c is close to a. 

V. RESULTS FOR HIGH AND MODERATE INTERFACIAL TENSIONS 

Typical snapshots of our nonequilibrium computer simulations are presented in Figure 
3 for 4 different times and for two different densities p* = pa 2 = 0.2, 0.4 and two differ- 
ent strengths of the driving force F* = Fa/ksT = 10,40. The parameter combinations 
investigated and the corresponding ratios of the line tension 7 and the driving force F are 
summarized in Table I. They are of the order 1 (for combination A-C) or a bit smaller (for 
combination D) such that the Rayleigh- Taylor instability wavelength is larger than or of the 
order of the mean interparticle spacing a, see again Table I for the ratio \ c /a. We call an 
interface tension "high" if \ c /a > 1 and "intermediate" if X c /a « I. This is in contrast to 
the case of low surface tensions (combination E) where \ c /a is significantly smaller than 1 
which is studied in chapter VI. 

The system was started with an equilibrated single interface situation. The initial inter- 
face is pretty smooth but carries capillary fluctuations on it. The density depletion shown in 
Figure 1, manifests itself as a clear void region in the interface reminiscent of a "forest-aisle" 
which is induced by the strong energy cost when two different particle species do meet. 

Then instantaneously the external driving forces were turned on forcing the particles to 
drift against each other. This first causes a local density increase close to the interface. 
Then interface undulations are getting more pronounced deforming the flat interface. This 
is the most efficient channel for different particle species to reverse their height. A typical 
characteristic undulation wave length can be identified from the simulation snapshot which 
we shall quantify later. Furthermore nonlinear effects such as interfacial overhangs ( "mush- 
rooms") can be seen (e.g. in Figure 3a for the largest time) until the mixture penetrates 
through eachother and the height reversal is complete. 

In order to quantify the structural signatures of the interface instability further, we 
consider the Fourier transform of the interface position as a function of time. In detail, let 
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h(x, t) be the interface position of a given configuration at a time t, and 



L 



h(k,t) — J h(x,t) exp (—ikx) dx (13) 
o 

the corresponding Fourier transform. We are interested in the time dependence of the 
averaged power spectral density which is defined as 

P(k,t) = (h(k,t)h*(k,t)) (14) 

where (...) now denotes an ensamble average over inital equilibrated configurations which 
dynamically evolved after a time t. If, at a given t, P(k,t) possesses a sharp maximum at 
k = k m , this implies that the interface will exhibit mainly undulations with a wave length 
of A m = 27c/k m . 

As an aside, let us remark that an alternative way of obtaining the line tension 7 is via 
the equilibrium capillary wave spectrum P(k, t = 0) which behaves as 

P(*.* = <» = ^ (15) 



for small wave vectors k 
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581 ] . We have checked that values for the line tension 7 as 
obtained from this formula are consistent with those obtained from Eq.(jHJ). 
We are now in a position to define a differential growth rate T(k,t) via 

r( *-') = ?(b)S p( *- ( ) < 16) 

A positive sign of T(k, t) implies that a mode of wave number k is growing at a time t while 
a negative sign means that - at a given time t - the mode is decreasing. 

In detail, the following numerical procedure was used to obtain the actual interface posi- 
tion h(x,t): For a given particle configuration at time t, we construct Voronoi cells around 
each particle and define the interface based on the associated Voronoi polygons. The poly- 
gon vertices that belong to cells of particles of both species A and B define a set of M* 
non-equidistant co-ordinates (h n , x n ) n=1 M ,. Now we eliminate the protuberance in the 
interface profile to get an functional interrelation for h(x). For further analysis we interpo- 
late (h n , x n ) n=1 M , by a cubic spline and eliminate sharp bends by using a simple low-pass 
filter 

h n ^h* n = (h n + h n+1 ) /2 (17) 
10 



where n — 1, . . . ,M and fiM+i = hi. This procedure results in a unique z-position h(x,t) 
of the interface at a given time t. 

Results for the differential growth rate T(k,t) are shown in Figure 4 as a contour plot in 
the plane spanned by the time t and the wave vector k. As in Figure 3, data for the four 
different parameters combinations A-D are presented. The white zero line summarizes 
points where T(k,t) vanishes and separates two regimes with growing (stable) modes and 
damped (unstable) wave numbers. Although this line is noisy, it can be read off from Figure 
4 that there is a whole band of wave numbers which are growing after an induction time of 
roughly tb- The lower limit of unstable wave numbers is significantly larger than the inverse 
simulation box length, while its upper limit is always smaller than a microscopic wave 
number of the order of 2ir/a. Between these two boundaries, roughly at the arithmetric 
mean of the two limits, there is a wave number with a maximal growth rate. !! The 
first basic observation is that for parameter combinations A-C where \ c /a is larger than 
1 the upper limit of unstable wave numbers is decreasing with time. !! This has to do 
with our equilibrated starting configuration. When the external field is turned on, there is 
a sedimentation-like process towards the interface which yields a local density increase at 
the interface. This behaviour can directly be seen in the simulation snapshots of Figure 3. 
Thereby the interface is getting stiffer as a fuction of time. We call this important effect self- 
healing of the interface, i.e. as a function of time the interface is getting less vulnerable with 
respect to short wavelength undulations. A possible destruction of the interface is efficiently 
blocked by a density accumulation. It is tempting to correlate this to an effective interface 
tension with a scaled density as obtained e.g. by the maximal total density p m close to the 
interface. This maximal density p m is a function of time. If one plugs this time-dependent 
density into the expression for 'y(p), one obtains a dynamically rescaled Rayleigh- Taylor 
expression (9). The corresponding wave number is also shown in Figure 4, see the open 
circles with error bars. Clearly the basic effect is encaptured by the dynamical rescaling as 
the classical Rayleigh- Taylor unstable wave number coincides well with the upper unstable 
limit. The agreement is in particular encouraging for parameter combinations A-C where 
\ c /a > 1. 

A special remark is in order for parameter combination D where A c is of the order of a. 
In fact, we estimate a microscopic wave number by ln^fp^ and data for this wave vector 
^\J~Pm are included in Figure 4 as crosses. Even in this case, the comparison between 
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the upper unstable wave number and the scaled Rayleigh- Taylor prediction is qualitative. 
This implies that the classic Rayleigh- Taylor criterion is astonishingly robust even close to 
molecular spacings. 

On the other hand, the maximal growth rate is seen as light region for small wave vectors. 
This wave number is slightly decreasing with time as well. We have further compared the 
wave number of maximal growth by calculating the second maximum of the height-height- 
correlation function which is defined as 



The position of the second maximum of C(k,t) is also included as a square- line in Figure 
4. It is increasing as a function of time towards larger wave numbers, i.e. smaller wave 
lengths. The increase is small and confined to a small "induction time" of the process. 
We attribute this to an initial process which has to do with the fact that due to the non- 
homogeneous capillary wave spectrum Eq. (jl5j) . undulations with larger wave lengths have 
initially a larger amplitude and have therefore more weight in C(k,t). This gives rise to 
the crossover to higher wave number after an induction time which is pretty sharp for the 
parameter combinations A and B. 

Finally we note that at higher k the growth rate is getting positive for larger times. This 
is due to the fact that the interface position is no longer sinoidal but starts to exhibit sharp 
parts close to overhangs. This is turn will generate higher Fourier modes to grow and this 
is what is indicated by the zero-line at higher k and larger times. 

We finally think that the zero-line at high k and smaller times (in Figure 4a and 4b) 
is statistical noise. Furthermore, a full comparison of the differential growth rates to a 
hydrodynamic approach requires a dynamically scaled viscosity and is left for future 
studies. 

VI. RESULTS FOR LOW INTERFACIAL TENSIONS 

Similar computer simulations were done as in the previous section but now in the regime 
where the line tension was low in the sense of the criterion of Eq. (jl2J) such that A c was 
significantly smaller than a. Simulation snapshots are shown for p* = 0.075 close to the 
estimated critical density in Figure 5. In this case, the interface is penetrable by the strongly 




(18) 
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driven particles and 7/i 71 = 0.002 and X c /a = 0.2, see the parameter combination E in Table 



re interface, the particles form lanes similar to the behaviour in additive 



I. 

In penetrating t 

mixtures in Refs. 25j, |26[. The width of the lanes is comparable to the correlation length 
of the mixture meaning that single worms of particles are formed (27 1. At the head of the 
lanes both spikes and extrusions are visible. A similar behaviour is found for even smaller 
densities below the critical density where the system is mixed in equilibrium. Hence the 
Rayleigh- Taylor criterion would predict a submolecular unstable wave length but the actual 
realized instability wavelength is a molecular correlation length of the system of the order 
of a. 



VII. CONCLUSIONS 



In conclusion, we have studied the onset of the Rayleigh- Taylor interface instability on 
length scales of interparticle distances for a compressible Brownian fluid mixture. Depending 
on the equilibrium interfacial tension, either a direct transition towards lane formation or 
the macroscopic Rayleigh- Taylor instability was found. An interesting dynamical effect is 
the self-healing mechanism for the interface which is produced by a density accumulation of 
particles driven against eachother which then causes an increase of the instability wavelength 
as a function of time. 

It would be interesting to study the dynamical process of phase separation by starting 
from a completely mixed system under the influence of the driving forces. Here one would 
expect a subtle interplay between phase separation kinetics which typically results in fractals 
and lane formation [59] . One should further investigate more extensively the dependence on 
the "sedimentation height" , i.e. how the data are affected by a larger initial D in ^-direction 
perpendicular to the interface. Moreover, it would be interesting to start the simulation by 
impressing a prescribed wave length as an interface undulation and checke its differential 
growth. A third set-up for another initial configuration which is typically used for the 
molecular dynmics simulations flO] would be to equilibrate first with respect to the external 
field with a fixed interface as induced e.g. by a thin hard platelet separating the two fluid 
phases. The effects of different starting configuration on the onset of the Rayleigh- Taylor 
instability will be left for future studies. 
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For simplicity, all our simulations were done in two spatial dimensions. In three dimen- 
sions, similar effects should persist. However, the characteristic wavelength has now two 
components parallel to the interface. The self-stabilizing effect of the interface and the two 
extreme limits should be similar as in two dimensions. The lane formation, for instance, has 
also be shown to be present in three spatial dimensions. 

Let us finally discuss some possible experimental verifications. Strongly non-additive 
mixtures which phase-separate into two fluid phases are found in colloid/polymer mixtures 
which exhibit colloid- rich and colloid-poor fluid phases for size ratio of about 0.5 or larger 
between the polymer and the cohoid 6QJ. Recently a phase-separating mixture of colloids 
and polymers was observed in sedimentation by Aarts and coworkers [17]. In particular, 
it was observed that large regions of phase-separated colloid-rich and colloid-poor phases 
exhibit a transition towards lane formation. The characteristic wavelength is of the order of 



that given by the Ray leigh- Taylor instability Eq. ffTTJ|) 61] . A full quantitative comparison 
has still to be performed. Another system were such lane formation has been seen and which 

n 

could be a candidate for a quantitative comparison is a xanthan-colloid mixture 62} . Finally 
we mention the fascinating possibility to study the mixture of complex plasmas involving 
dust grains for interface instability in real space 63, 
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FIGURE AND TABLE CAPTIONS 



Figure 1: Averaged density profiles pa{z) and pb(z) of the two particle species and 
anisotropy Pjy — Pt of the pressure tensor as a function of the distance z perpendicular to 
the interface for A = 1.6 and for four different densities: pa 2 = 0.4 (solid lines), pa 2 = 0.3 
(dashed lines), pa 2 = 0.2 (dot-dashed line) and pa 2 = 0.15 (short-dashed lines). 

Figure 2: Reduced line tension jk B T as a function of reduced bulk density pa 2 for two 
different nonadditivities A = 0.8 (solid line) and A = 1.6 (dot-dashed line), (a) Away from 
the critical point, (b) Closer to the critical point. The estimate for the critical density are 
shown as arrows. 

Figure 3: Typical simulation snapshot for an interface in nonequilibrium with an external 
drive. There is a sharp interface between an A-rich and a B-rich fluid phase. The starting 
configuration was an equilibrated interface at z — 0. Four different times are shown as given 
in the legend, a) pa 2 = 0.2 and F* = Fa/k B T = 10, b) pa 2 = 0.2 and F* = Fa/k B T = 40, 
c) pa 2 = 0.4 and F* = Fa/k B T = 10, d) pa 2 = 0.4 and F* = Fa/k B T = 40. 

Figure 4: Contour plots of the growth rate T(k,t). Zero growth is shown by the white 
line. The rescaled classical Rayleigh- Taylor wave length 2ir/\ c is shown by the circles. Here 
is time-dependent maximum p m of the total density at the interface is used. The statistical 
error stems from the uncertaincy of this density. The position of the second maximum in 
the height-height correlation function are the squares. The corresponding microscopic wave 

1/2 

number as defined by lixpm and gives an estimate for the threshold to where a microscopic 
length scale are realized. 

Figure 5: Typical simulation snapshot for an interface in nonequilibrium with an external 
drive. There is a sharp interface between an A-rich and a S-rich fluid phase. The starting 
configuration was an equilibrated interface at z — 0. Four different times are shown as given 
in the legend. The parameters are pa 2 = 0.075 and F* = Fojk B T = 80. 

Table 1: Summary of the 5 parameters combinations studied in the paper. The corre- 
sponding ratios j/F and \ c /a are also given. 
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